Dynamic structure selection and instabilities of driven Josephson lattice in 

high-temperature superconductors 
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We investigate the dynamics of the Josephson vortex lattice in layered high-T c superconductors at 
high magnetic fields. Starting from coupled equations for superconducting phases and magnetic field 
we derive equations for the relative displacements [phase shifts] between the planar Josephson arrays 
in the layers. These equations reveal two families of steady-state solutions: lattices with constant 
phase shifts between neighboring layers, starting from zero for a rectangular configuration to tt for a 
triangular configuration, and double-periodic lattices. We find that the excess Josephson current is 
resonantly enhanced when the Josephson frequency matches the frequency of the plasma mode at the 
wave vector selected by the lattice structure. The regular lattices exhibit several kinds of instabilities. 
We find stability regions of the moving lattice in the plane lattice structure - Josephson frequency. 
A specific lattice structure at given velocity is selected uniquely by boundary conditions, which are 
determined by the reflection properties of electromagnetic waves generated by the moving lattice. 
With increase of velocity the moving configuration experiences several qualitative transformations. 
At small velocities the regular lattice is stable and the phase shift between neighboring layers 
smoothly decreases with increase of velocity, starting from n for a static lattice. At the critical 
velocity the lattice becomes unstable. At even higher velocity a regular lattice is restored again 
with the phase shift smaller than tt/2. With increase of velocity, the structure evolves towards a 
rectangular configuration. 
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I. INTRODUCTION 

The intrinsic Josephson effectEl in strongly anisotropic high-temperature superconductors, such as Bi2Sr2CaCu20 :E 
(BSCCO), has been a subject of intense experimental study (see review articlesnB). Stack of low-dissipative junctions 
formed by atomic layers represents a nonlinear system with unique dynamic properties. In particular, this system 
possesses a rich spectrum of the electromagnetic plasma waves. These waves can be excited via the Josephson 
effect, which makes these materials convenient voltage-to- frequency converters and gives strong potential for practical 
applications. 

A magnetic field applied along the layers forms the lattice of Josephson vortices. Transport properties of this 
state are determined by the dynamics of the. Josephson vortex lattice. Coherent motion of Josephson vortices in 
BSCCO has been observed experimentallyotl The motion of the Josephson lattice generates alternating electric 



ne can expect a strong resonance emission 
. In a voltage-biased transport experiment 
I, while in a current-biased experiment the 



fields and currents which are coupled to electromagnetic plasma waves, 
when the velocity of the lattice matches the velocity of the plasma wavi 
this should be seen as a resonant enhancement of the transport current! 
resonance should be seen as a voltage jump. In conventional long Josephson junctions this effect is known as Eck peak 
(Eck step)li3 (see also recent article O and references therein) . However, in layered superconductors this effect has 
qualitatively new features. The velocity of a plasma wave propagating along the layers strongly depends on the wave 
vector along the c direction. On the other hand, the traveling Josephson lattice is coupled only to the plasma modes 
with the wave vectors given by the reciprocal wave vectors of the lattice, i.e., the^mve vector of the resonant mode is 
determined by the structure of the Josephson latticecn. In early theoretical workErta it was assumed that the moving 
lattice preserves its static structure (stretched triangular lattice) and the resonant velocity of the lattice is determined 
by this structure. We will show that the lattice structure experiences a nontrivial evolution with increase of velocity. 
Moreover, a regular lattice state always becomes unstable at a certaiEL-iplocity. Therefore, finding conditions for the 
resonance emission is a nontrivial problem. Recent numerical studiestJiiJ revealed a very rich dynamic evolution of 
moving Josephson structures, which includes regular as well as chaotic solutions. 

In this paper we analyze the dynamics of the Josephson lattice at strong fields, in the regime of strongly overlapping 
vortices. In this regime Josephson coupling can be treated perturbatively. This strongly facilitates the analytical 
analysis of the dynamic structures. Major results have been already reported in the Letter |l5|. A similar approach has 
been, used— to calculate resonant current-voltage dependencies in conventional long Josephson junctions in magnetic 
field .IIEIIIj'EJ This approach also describes Josephson dynamics at arbitrary magnetic field when all junction are driven 
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into resistive state by transport current. We will focus on the dynamics in the case of large size samples and will not 
consider resonances related to finite-size effect (Fiske resonances) . The classification of Fiske respaant frequencies for 
layered superconductors in the case of a moving rectangular, lattice has been done by Kleiner Jl_ The experimental 
observation of these modes has been reported by Irie et. al.E_l ______ 

The phase dynamics in layered superconductors is described by the coupled sine-G-D-an-c quationsM___£l__ Two 
equivalent approaches have been used to derive these equations. In the first approacl_li_'L_ the layered superconductor 
is modeled as a multilayered S-I-S-I-. . . system and the parameters of the equations are expressed via the geometrical 
parameters and the bulk parameters of the superconducting material from which the superconducting layers are 
prepared. The second approacl_Jtr_, more natural for atomically layered superconductors, is based on the Lawrence- 
Doniach model. In this approach the layers are treated as two-dimensional entities from the very beginning. We will 
use the second approach. 

The dynamic phase diagram strongly depends on dissipation mechanisms. Moving Josephson vortices generate both 
in-plane and out-of-plane electric fields which induce dissipative quasiparticle currents. Usually, only dissipation due 
to the tunneling of auasiparticles between the layers is taken into account in consideration of the dynamics of the 
Josephson vortices. c3tB However, in high-T c superconductors the in-plane component of the quasiparticle conductivity 
a a b is strongly enhanced in the superconducting state as compared to the normal conductivity due to the—reduction 
of the phase space for scattering. This enhancement was observed in YBCI2CU3O713 and «i-SS'CCC_l. The c- 
axis component, er c , monotonically decreases with temperature in the superconducting state.l_iE_ The anisotropy of 
dissipation u a bl^c becomes larger than the superconducting anisotropy 7 2 (here 7 is the anisotropy of the London 
penetration depth) . This leads to the dominance of the in-plane dissipation in the dynamics of the Josephson lattice 
in a wide field range. _] 

Starting from the the coupled sine-Gordon equations, we derive equations for the relative displacements between the 
moving Josephson planar arrays in the layers (or, equivalently, for the phase shifts between the Josephson oscillations 
in the different layers). This allows us to reduce the original two-dimensional problem to an one-dimensional one. 
Dynamic equations reveal two families of steady-state solutions: lattices with constant phase shifts between neighbor- 
ing layers, starting from zero for a rectangular configurations to 7r for a triangular configuration, and double-periodic 
lattices with the average phase shift 7r/2. We analyze the stability of the regular lattices and show that they exhibit 
several kinds of instabilities including the long- wave length shear instability and the short-wave length instability with 
respect to alternating deformations. A specific lattice structure at given velocity is selected uniquely by boundary 
conditions. The numerical investigation shows that at small velocities the lattice experiences a smooth evolution of 
the structure with the phase shift between neighboring layers decreasing from -n at zero velocity to smaller values. 
At a certain velocity the lattice becomes unstable and the system switches into an oscillating or chaotic state. At 
even higher velocities the steady-state regular lattice is restored again. The phase shift for this rapidly moving lattice 
is smaller than 7r/2 and continues to decrease with increasing velocity, i.e., the structure evolves smoothly towards 
a rectangular configurations. The current-voltage characteristic is non-monotonic in the region of the instabilities, 
which means that two different coherent steady-states coexist within a certain current range. Within this current 
range the system can also be in the phase separated states, in which it is split into two (or more) domains moving 
with different velocities separated by a phase defect. 

The paper is organized as follows. In Section || we derive the coupled dynamic equations for the in-plane fields 
and phase differences. In Section III we derive the dynamic equations for the relative displacements [phase shifts] 
between the Josephson planar arrays in layers. In Section |y] we study the coherent steady-states [regular lattices] 
and calculate the excess Josephson currents and Poynting vectors for such states. In Section [v] we investigate the 
stability of the steady-states. We consider in detail two important particular cases of instabilities: the long-wave and 
short-wave instability. We also compute numerically the stability phase diagrams in the plane structure- Josephson 
frequency. In Section VI we consider boundary conditions for top and bottom boundaries, discuss the mechanism of 
structure selection by the bou ndar ies, and numerically calculate the evolution of the lattice structure with increase 
of velocity. Finally, in Section VII we discuss the mechanism of multiple branches in the current- voltage dependence 
due to the dynamic phase separation: spontaneous splitting of the system into the rapidly and slowly moving regions. 



II. DYNAMIC EQUATIONS 

The dynamics of the moving Josephson lattice can be described by coupled equations for the phase differences and 
magnetic field. These equations can be derived from Maxwell's equations expressing fields and currents in terms of 
the gauge invariant phase difference between the layers 9 n = <fi n +i — 4> n — ^^-A z and the in-plane superconducting 
momentum p n = V x cj) n — t^A x (see, e.g., Ref. |9|). Consider a layered superconductor in a magnetic field applied 
along the layers (y -direction) with transport current flowing along the c-axis (z-dircction). A local magnetic field H n 
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between the layers n and n + 1 can be expressed as 



H n(x) = ^7 ( ?7 - Pn+l +Pn)- ' 1 1 



27rs \ <9o; 

The components of electric held can be approximately represented as 

~ 2nc dt ' z ~ 2^ CS 9* • 1 j 

The components of electric current, and j z , include the quasiparticle and superconducting contributions 

Jx_,Tafc 27rc dt S^A 2 /"' [ ' 

' jjsm6» n . (4) 



27TCS dt 

where a a b and tr c are the components of the quasiparticle conductivity, jj — 8 ^" A2 is the Josephson current density, 
and A a b and A c are the components of the London penetration depth. Using the above relations we rewrite the z- and 
i-component of Maxwell's equation -^j + ^ = V x H as 

2a c $ d9 n 4tt . e c $ d 2 #„ 9# n 

flf + T JjSm0n + 2^^ = -fe' (5) 

2cr aD *o dp n $0 #™ - #n-l ,„s 

2 «T + o \2 P» = • ( 6 ) 

c 2 dt 27rA ab s 

In the second equation we replaced dH/dz by discrete derivative (H n — H n _i)/s and neglected the in-plane displace- 
ment current (typical frequencies are assumed to be much smaller than the in-plane plasma frequency c/A ao ). Taking 
discrete derivative of the second equation and using relation (|l|), we rewrite the second equation in terms of 9 n and 

f 4%a ab d , 1 ^ f$ 39 n v \ H n+1 +H n _ 1 -2H n 

V c 2 dt A 2 J {2tts dx n ) - s 2 ■ [ ' 

Equations (||) and (^) describe the dynamics in terms of phases 6 n (x,t) and field H n {x,t). They are equivalent to 
equations derived in Refs. ||j9[ In the case of negligible in-plane dissipation (cr ao = 0) these equations can be reduced 
to coupled system of equations for the phases 9 n , which were used in Refs. |i|],[ij|20]j2l]. Equations (|) and (0) provide 
a simplified description of the phase dynamics. In real systems this dynamics can also be influenced by the charging 
effectsEj and the quasiparticle imbalances. However, at present there is no experimental information, which would 
allow to describe these effects quantitatively. 

To simplify the analysis of dynamic equations (^|) and (Q), we introduce the reduced coordinate u = x/js, time 
t = ujpt, and magnetic field h n = H n 2-K^\ 2 /$o- We also introduce the reduced penetration depth I = ^ut and the 

dissipation parameters (reduced conductivities)^ v c = |^pS v a b — i7r(Ta ^ " p , which are related as v c jv ab = a c j 2 /a a b- 
In terms of these variables Eqs. (||) and (Q) acquire a much simpler form 

a2f) d6 n . „ dh n n 

v c ^ + sm8 n -—^=0, (8) 



dr 2 dr du 



2 1 \ , d9 n : d ( 39 n h n 



^-p) h « + l£ + V ° b d-r{l£-f) = °< (9) 

which we will use in our analysis. Here V 2 denotes the discrete Laplace operator, V 2 /i„ = h n+ i + /i n _i — 2h n . 

In certain special cases the quasiparticle conductivity can be directly related with the Josephson current.. The most 
known-example of such relation is the Ambegaokar-Baratoff formula for conventional Josephson junctions.Ej Latyshev 
et. alEB have found a similar relation for BSCCO at small temperatures, jj = ira c A / (es), with A being the 
maximum gap. This relation was derived assuming d-wave symmetry of the order parameter and coherent interlayer 
tunneling. Using this relation we obtain a simple estimate for the c-axis dissipation parameter v c at small T, 

v c « -^-,at T -> 0. 
2ttAo 
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Taking parameters, typical for slightly underdoped BSCCO at T w 50K: 7 = 500, \ a b = 240 nm, s = 15 A, p a b = 
cr^ 1 = 50 ^Sl-cm, and p c = a^ 1 = 600 il-cm, we obtain estimates ^ a h « 0.1 and v c w 0.002, which we will use in the 
numerical computations. As we can see, typically the in-plane dissipation is much stronger than the c-axis dissipation. 
Models, which take into account only the c-axis dissipation channel, are not suitable for the quantitative comparison 
with experiments in BSCCO. Meanings, definitions, and practical formulas for the dimensionless parameters are 
summarized in Table 1. 



III. EQUATIONS FOR RELATIVE DISPLACEMENTS [PHASE SHIFTS] BETWEEN THE JOSEPHSON 

PLANAR ARRAYS. 

At high magnetic fields H > $0/(5. 57s 2 ) (h > Z 2 /5.5) strongly overlapping Josephson vortices fill all the layers.El 
In this regime the Josephson coupling can be treated perturbatively. Without the Josephson coupling (i.e., dropping 
sin#„ in Eq. (g)) the phases 9„ change in space and time as 

0<W (r, u) = luet + k H u + <t> n (10) 

where the frequency uje is set by the electric field E z , loe = 2izcsE z / (<f>aL) p ) and the smooth function kn{u) is 
connected with the local magnetic field as d(knu) / du = h/l 2 . In general, the local field h varies in space due to 
self- field of the transport current. In the following we consider a situation when the self-field of current is much 
smaller then the applied field and neglect its contribution to h. In this case kniu) can be approximated by a 
constant, fejj = h/l 2 = 2nH~/s 2 /<&o. Therefore Eq. ( [l0| ) represents the phase modulation moving with the velocity 
ve = uJE/ku = cE z /H. The phase shift <f> n characterizes a relative coordinate of the Josephson lattice in the layer 
n. Change of the phase shift Stfi n corresponds to the displacement of the lattice at distance S(p n /kjj- Without the 
Josephson interaction the system would be degenerate with respect to the phase shifts <j> n . The Josephson coupling 
eliminates this degeneracy. It either leads to the slow dynamics of the phase shifts </>„ or selects a specific steady- 
state structure. We will derive a closed set of equations governing the dynamics of the phase shifts <^„(r). This 
can be done using an expansion with respect to the Josephson currents and averaging out rapidly oscillating degrees 
of freedom. Since the compression deformations of the lattice have much higher stiffness in comparison with the 
shear deformations, one can neglect the explicit coordinate dependence of </>„. In the first order with respect to the 
Josephson current, the phases Q n (r,u) and field h n (r,u) acquire the oscillating corrections, 9 n (r,u) and h n (r,u), so 
that they can be represented as 

#„(t, u) = luet + k H u + 4> n (T) + 9 n (r,u), (11) 
h n (r,u) = h(u) + h n (r,u), (12) 

with (d(f> n /dT) n — 0. Substituting 9 n (T,u) into Eq. (||) and averaging it with respect to u we obtain the current 
conservation condition 

Here ij n = (sin0 n ) u is the local Josephson current between the layers n and n + 1 averaged with respect to the 
in-plane coordinate u and ij — (sm9 n ) = ((sin is the Josephson current averaged over the layer index. The 
last current obeys the macroscopic equation v c loe + ij = dh/du. In the following, we will assume that d^njdr <C loe- 
To obtain a closed set of equations we have to relate the local currents ij n with the phase shifts <fi m and their time 
derivatives d<p m /dT. For this we have to find the oscillating phases and field and average sm9 n with respect to the 
rapid space variations. The total Josephson current ij n = ij n [<f> m ,d<p m /dr] consists of the nondissipative (i' Jn ) and 
dissipative (i" n ) components 

*Jn=ij n + ij n - (14) 

The nondissipative component i' Jn determines interaction between the phase shifts <p n . To the lowest order with 
respect to d4>n/dr it depends only on the instantaneous configurations of the phase shifts i' Jn = i'j n [4>m]- The 
dissipative component i'j n gives an additional contribution to the dissipation coming from the Josephson coupling, 
and it is proportional to d<p n /dT for slow time variations. The first term is determined by the first order oscillating 
corrections to the phases 9 n {u) and fields h n which can be found from equations 
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d 2 9 n ae n dh n 

-g-^- - u c~g~T + ~g~~ = sm \ w et + k H u + <pn) , (15a) 

v «~^j^ + ^7 + ^U7"^) =0 - (15b) 



The Josephson term in the first equation acts like a source generating the electromagnetic wave with the frequency loe 
and in- plane wave vector kn- To obtain solution of this equation we introduce the phase and field response functions, 
G(n,m) = G(n,m; k 7 ui) and H(n,m) = H(n, m; fc, w), defined as solutions of equations 

(lj 2 — iv c uS) G(n, m) + ikH(n, m) = (5„ !m , (16a) 
ik (1 + iv ab uj) G(n, m) + (v 2 n - 1 + ™ abLJ ) #( n , m ) = . (16b) 



i 2 

These complex functions give the oscillating phase and field in the m-th layer generated by the oscillating cur rent in 



the n-th layer. For bulk response both functions depend only on the difference n — m and the solution of Eqs. (16a, b) 
can be obtained by Fourier transformation. For the function G(n — m) this gives 

/dq 
— exp(iq(n-m))G(q,k,ui), (17) 
27T 

2 k 2 (1 + iv ab uj) 



Q(q, k, u>) = I u> — iv c ui 



2(1 - cosg) + (1 + iv ab u>) /l< 



This function will play a very important role in the following derivations. In general, the parameters v c and v ab in 
the last formula may depend on frequency due to the frequency dependence of the quasiparticle conductivities. The 
pole of the function Q(q,k,u) gives the spectrum of the plasma waves and their damping. The spatial dependence 
of the function G(n — m; k, u>) is determined by the complex wave vector q + = q+(uj, k), which describes propagation 
and decay of the plasma wave along the z axis for the given frequency lo and in-plane wave vector k, 



Sn k 2 (1 + iv ab uj) exp(zg + |n| ) 

iv c u> (lu 2 - iv c uj) 2 2isinq + 



G(n; k,u) = : , 2 \ ^ ab 2 J 2 „ ■ ( 18 ) 



k 2 (l + w ab ui) l + iv ab w 

cosq + = l-—— ; - + — , with Im g+ > 0. (19) 

2 (lo 1 — iv c lo) 2l z 



(20) 



Now we can represent the solution of Eqs. (15a,b) for the oscillating phase 9 n (T, u) in the form 
n (r, u) = Im G(n — to; kn,0JE) exp (iloet + iknu + i<j> m ) 

m 

and calculate the nondissipative current i' Jn [<j> n ] in Eq. (|l^) : 

i'jni&n] = (0„(t,u) cos (ll> e t + k H u + 4> n )) u (21) 
= 2 Im t G ( n _ m; W£ ) exp ^™ ~ ^ m ^] ■ 

m 

This term determines the nonlocal dynamic interactions between different (f> n . As follows from Eq. (|l8|), the range 
of interaction is given by the length L s (ku ,we) = (Im[(7-|_(fc#, wg)]) -1 . At = interactions extend only over the 
nearest neighbors. The range of interactions increases with increase of lattice velocity. When the velocity exceeds the 
minimum velocity of propagating electromagnetic wave, the range of interactions is limited only by dissipation. 

To describe the slow phase dynamics, one has to find also the dissipative component of the Josephson current 
ij n [(j> n ,d<f> n /dT]. It is determined by the oscillating contributions to the phases and fields of the order of n d<f) n /dT 
which we notate as Qdn and hdn- Representing 9 n and h n in the complex form 6 n , h n oc exp (iujet + iknu), we obtain 
the following equations for 9d n and hdn 

(lo 2 e - iv c oj E ) dn + ik H h dn = -TT~ (-2w E + iv c ) &n, 

OT 

Y72u \ + W ab LO E . -7 / -i . , a d(j> n ( . K i n\ 

V n h d n J2 hdn + ^ + lv abUE)Vdn = I ™ ab — + V ab k H 8 n I . 
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Solution for 6d n can be represented as 



'dn 



F (n — m,ni — to) ■ 



<9t 



• exp {iuo E T + ik H u + icp m ) 



where the kernel F (n, to) = F (n, m; Uh^e) is denned by relations 

F(n,ni) = f y^^J 7 (q,qi)exp(iqn + iqini) , 
J (27r) 



-2wi 



2(1 — cosgi) 



2(1 — cosg) 



2(1 — cosgi) 



I 2 



G(q)G(qi 



(22) 
(23) 



with G{q) = G(q,k,uj) being the response function defined by Eq. (|l"7j). Using this solution we can now obtain the 
term i"„[0„, #„/dr] 



i"„[0„, d^n/dr] = (6> dn (r, u) cos (w B r + + </>„))„ 



(24) 



_F (n — m, n\ — to) 



Ot 



exp (i (0 m - 4> n )) 



Substituting averages (HI) and (||) into Eq. (|l|) we finally obtain: 

^n-m ^p 1 + ^ ^ Im [G(n - to) exp ( 



071 ))] = *J, 



(25) 



where the nonlocal viscosity matrix v n - m [4>j\ depends on the instantaneous phase configuration and is given by 

v n - m [4>j] = v c 8 n - m + ^ J~] Im [F (n - m, m - j) exp (i (<f>j - </>„))] (26) 



and the average Josephson current ij is given by 

ij = - y (Im [G(n - to) exp (-i (</>„ - <f> m ))]) T 



(27) 



These equations describe the slow dynamics of the phase shifts <j) n imposed by the Josephson coupling. Our further 
analysis is based on these equations. Eqs. (|25|) and ( p6f ) are valid far away from the top and bottom. This means 
that in a stack consisting of a finite number of junctions N, these equation have a region of validity only if N much 
larger than the decay length of electromagnetic wave along the c direction L s (kH ,ue) — (Im[g + (fcf-f , uje)])" 1 , which 
determines the range of the phase interactions. In Fig. g we plot the dependence of the decay length L s on the 
Josephson frequency loe for several values of ku using parameters v a b = 0.1 and v c — 0.002. For small fields one can 
see a significant increase of L s when loe exceeds the minimum frequency, Lo m i n = kn/2, for the electromagnetic wave 
with the in-plane wave vector ku- Nevertheless for stacks, containing 50-100 junctions, the bulk behavior is expected 
in a wide range of fields. The same stack can be either in the bulk regime or in the finite-size regime depending on 
electric and magnetic fields. 



IV. COHERENT STEADY-STATES AND RESONANT I-V DEPENDENCIES 

In steady-states we have d4> n /dr = and the static phase shifts <j) n far away from the boundaries obey equations 

i ^2 Im [G{n - to) exp (-i (0„ - (f> m ))} = ij. (28) 

rn 

These equations have a trivial solution 



4>n = Kn, 



(29) 



G 



corresponding to a regular lattice. In this solution the planar arrays in the neighboring layers are shifted by the 
fraction k/2it of the lattice constant (see Fig. 0a). In particular, n = corresponds to a rectangular lattice, and n = it 
corresponds to a triangular lattice giving a ground state at oje = 0. The oscillating phase (f20|), 8 n = 9 n (K, kn ,o->e), 
and the average Josephson current (|27|), ij = ij(K, kH,^E)i f° r the moving regular lattice are given by 

6 n = Im [Q(k, k H ,u E ) exp i (u b t + k H u + nn)] , (30) 
ij = -Im[g(K, k H ,uj E )} . (31) 

In the case of relatively weak in-plane dissipation, v a b^E "C 2Z 2 (1 — cos k) + 1, the Josephson current acquires a simple 
resonant dependence on the Josephson frequency oje (electric field) 

■ _ 1 uev{k) ,„ 9 X 

^(^- W 2( K )j +( W£;! y( K )) 

where 



oj p (k) = 

is the plasma frequency at the wave vector k and 

v (k) = v c 



k 



H 



^(l-COS^ + l// 2 

2 (1 — cos k) k 2 H v a i> 
(2 (1 - cosk) + ±) 2 



is the K-dependent dissipation parameter. When the Josephson frequency oje matches the frequency of the corre- 
sponding plasma wave lo p (k), a resonance enhancement of the current is expected. Using this result we can represent 
the current-voltage characteristic as 



1 E$v(k)/v c 

2 {El - E?( K )) 2 + {E p E z v{n)f 



3z = °cE z [ 1 + - — — P J '' - - o , (33) 



where 

E r {n) 



Hs/X ab 



^T c ^/2(1- cos K ) + s yXi b 



is the resonance electric field and E p = $>oU> p / (2ncs) is the electric field at which-the Josephson frequency matches 
the plasma frequency. This expression is similar to Eck peak for a single junction.!!!] However, in contrast to a single 
junction, both the resonance frequency and the dissipation parameter depend on lattice structure. An important 
scale of the electric field is set by the minimum frequency of the propagating plasma modes, corresponding to the 
triangular lattice 



E min = E t (-k) 



Hs 



4A„ 



For BSCCO [\ a b = 220 nm, and s = 15 A, e c = 11) this field typically corresponds to the voltage drop of rj 0.48 mV 
per junction at H = 1 tcsla. At fixed electric field the triangular configuration (k — n) gives the highest dissipation 
and the highest current, while the rectangular configuration (k = 0) gives the lowest dissipation and the smallest 
current. There are two reasons for this property. Firstly, the triangular configuration has the smallest resonance 
frequency oj p (tt) « fcfj/2 and the rectangular configuration has the largest resonance frequency uj p (0) = Ikn exceeding 
oj p (t:) by the factor 21 = 2A a t/s ~ 300. Secondly, the in-plane contribution to dissipation is completely absent for the 
rectangular lattice and it is maximum for the triangular lattice. Eq. ( |33| ) can be used for the direct comparison with 
experiment only at small electric fields, where the structure is close to triangular one. At higher electric fields the 
lattice structure (parameter k) changes with the electric field. This evolution of the lattice structure is determined 
by the boundary conditions and will be considered below. 

The motion of the regular Josephson lattice generates a traveling electromagnetic wave inside the superconductor. 
The energy flux in this wave is given by the Poynting vector 
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P = ^-[ExH] 



where E and H are the oscillating components of the electric and magnetic fields. Using expression for the y component 
of the magnetic field, H y , and for the components of the electric field, E x and E z , 



<j> 



27TCS7 
I'KCS 



ik H (l + iv ab LU E )g(n, k H ,LO E ) 
2(1 - cosk) + (1 + iiy a bLOE)/l 2 
k H uj E (1 - exp (-£«;)) £?(«;, k H ,uj E ) 



2(1 - cos k) + (1 + w ab uj E )/l 2 
Im [ui)eG{k, kn , u>e) exp (i (fcjjcc + /tn + u E t))} 



exp (i (fcjja; + cn + uJ E t)) 
exp (i (feiya; + kii + we*)) 



we obtain the components of the Poynting vector, 

o; g fc g (2(1 - cosk) + (1 - v 2 ab 0J 2 E ) / 1 2 ) 
\2{l-cos K ) + {l + iis ab LU E )/P\ 2 



P. 



p, 



uj E k 2 H (sinK — v ab uj E (1 — cos ft)) 
|2(l-cosK) + (l+w ab w jE )// 2 | 2 



\G(k, k H ,oj E )\ 2 , 



\G(K,k H ,W E )\ , 



with 



Pab 



32tt 3 A 2 s7' 



Pc = Pab/l- 



(34) 
(35) 

(36) 



For the typical parameters of BSCCO (A = 200nm, 7 =500, aj p /2ir — 150 GHz), the scale of the in-plane Poynting 
vector can be estimated as, P ab w 125 W/cm 2 . 

At small electric fields the lattice is triangular and Eq. ( p3| ) gives the linear flux-flow conductivity. We represent it 
in the form convenient for comparison with experiment 



ff ~ 
ai J ~ Or 



1 / $0 



2 \ttH 1 s 2 




&ab ( ^0 



27 2 \irH"fs 2 



(37) 



We neglected here small terms ~ I /I 2 . This result, valid at H > Qo/ttjs 2 , was first derived in Ref. ||. An impor- 
tant feature of this dependence is the H~ 2 term coming from the in-plane dissipation. In high-T c superconductors 
suppression of scattering leads to a strong enhancement of Jthe_in-plane quasiparticle conductivity, which is seen 
as a pronounced peak in the temperature dependence of (TafcEja. On the other hand, a c monotonically decreases 
when the temperature goes down. The inequality a ab 3> er c 7 2 holds almost everywhere in the superconducting state. 
In this case Eq. (^7|) predicts the quadratic dependence of the flux flow resistivity pfj cx H 2 in a wide field range 
<f>o/7T7s 2 < H < ■J <r ab /<Tc7 2< i > o / 7i"7s 2 £3 In contrast, for the case of dominating c-axis dissipation, the flux-flow conduc- 
tivity at high fields coincides with the c-axis quasiparticle conductivity and it is field independent. Eq. ( |37j ) describes 
very well the field dependence of the flux-flow resistivity in BSCCO. 

Another important case, which is realized at high currents, is the rectangular lattice with k = 0. In this limit Eq. 
(E53|) significantly simplifies 



j z = a c E z 



1+i 



2 (E 2 - E 2 r (0)) 2 + (E p E z v c ) 2 



(38) 



where E r (0) = Hj y/eZ is the largest resonance field. Note that in this limit the in-plane dissipation does not influence 
any more the lattice dynamics. A pronounced resonant feature at E z — E r (0) can be observed only for the case of 
very small c-axis dissipation. The peak amplitude is comparable with the quasiparticle current a c E z if 

H < <S> /(2nv c \ c s). 

For typical BSCCO parameters this corresponds to fields < IT. Another obvious condition for observation of the 
resonance is that the voltage drop between the neighboring layers at the resonance field sE r (0) must be smaller than 
the gap voltage A /e, which gives 
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H < Ha= V 

es 

Using Ao = 25mV and e c = 12 we obtain Ha ~ 2kG. 

In addition to the periodic lattices, Eq. (|2^) allows for the double-periodic solutions of the form 

4>n = |rc + u(-l)" 

This can be shown by the direct substitution of this expression into Eq. (|2|) using the following identities 
exp {-iu ((-!)" - (-ID) = 1 + (1 ~ ( ~ ir " m) (exp (-i(-l)"2u) - 1) 



^2 G ( m ) ex P (~*^ m ) = ^2 G ( TO ) exp y\ 

m 

ce is sketched 

parameter u, ij(kn,u}E) = \ Im , k 



to 

2 



Such double-periodic lattice is sketched in Fig.JIlb. The current ij for these states does not depend on the modulation 



V. STABILITY OF COHERENT STEADY-STATES 



A. Equations for small deviations from regular lattice and stability criterion 

To investigate stability of the moving regular lattice we consider a solution of the form 4> n {T) = un + u n (r) and 
obtain equations for the small deformations u n (r) <C 1 

vS^r~ + \ ^ Im \ F ( n - n ii K ) ex P (wt(»i - n))] (39) 
or 2 ' or 

m 



— Rc [G(n — to; uoe) exp [in(ra — n))] (u m — u n ) = 0. 



2 

m 

Looking for solutions in the form of plane waves 

u n = Re[w g exp(a(g)r + iqn)], 
we obtain for the eigenvalue a(q) = a(q, K, kn, Wb) 

1 g(« + g)+g*( K -g)-2Re[g( K )] 

a W = ~7 1 i \t( _i_ ^ TT( I YT ( 40 ) 

4 + 47 (k + 9, k) - ^* (-« + q, K)\ 

where the functions G(q) and T(q, qi) are given by Eqs. (|rj| ) and (p3j). The lattice is stable if there is no exponentially 
growing solution in the whole g-interval < q < it, i.e., 

Re[a(q)] < 0, for < q < it. 

The instability is characterized by the wave vector qi, at which Re[a(g)] becomes positive for the first time. We 
consider in detail the two important special cases: the long- wave instability at qi — and the instability with respect 
to alternating deformations at qi = it. 

According to the generally accepted classification (see, e. g., Ref. HJ), one can distinguish two kinds of instability: 
an absolute instability, for which initial perturbation exponentially grows at any point and a convective instability, 
for which growing perturbation is drifted away so that it decays with time at fixed point. 



B. Long- wave length stability 

Formally, the condition of the long-wave stability can be obtained from Eq. (|o|) at q — > 0. However, for better 
understanding of the instability mechanism, we will obtain it directly from the equation of motion for the elastic 
deformation of the lattice (|39|). To derive this equation, we add to the solution (g9|) slowly changing in space and 
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time function u(r,n), <j) n — nn + u(r, n), and replace the discrete index n by a continuous variable z. Substituting 
this expression into Eq. (^B|), expanding it with respect to u(r, z), and performing a gradient expansion, we obtain 
equation for the elastic deformation u(t, z): 



du 
"IT +C 



with 



dr dzdr 



a% = ax (k, k Hl uj E ) = -z Im 



du d 2 u 

a X ~7T~ + a 2TT^J 

OZ OZ A 



dG(K,k H ,uj E ) 



= 



a 2 = a 2 (n,k Hl uj E ) 



Re 



8 k 



d 2 g{n,k H ,uj E ) 



v = v(k, k H , uj e ) — v c + - Im 



8k 2 
dQ(n,k H ,uj E ) 



C = C(K,k H ,u E ) 



Re 



8oje 



di(n,k H ,ui E ) 
duj E 



dq 



(41) 

(42a) 
(42b) 
(42c) 
(42d) 



9=0 



In the derivation of the formula for v (42c) we used the relation J-(n, K,k,<J) = dQ(K,k,Lo)/du>. Note that the 
coefficients ax and C vanish for the symmetric lattices with k = 0, tt. Substituting a solution in the form of a plane 
wave u(r, n) cx exp (a(g)r + iqn) into Eq. (41) we obtain for small q 



a(<l) ~ ~ (-iaiQ + («2 - aiC/f) g 2 ) • 
This gives the stability condition Re[a(g)] < 

0-2 — dxC/ v < 0; v > 0- 



(43) 



(44) 



The first of these conditions means a positive shear stiffness, while the second condition is equivalent to the condition 
of the monotonic I-V dependence at fixed k, di/dui E > 0. On the other hand, formally, the lattice is also stable when 
both inequalities are opposite 



(45) 



The case ax 7^ corresponds to a convective instability. This means that only for symmetric lattices the long-wave 
instability is an absolute instability. Below we consider several important special cases of the long-range instability. 

For the important particular case of the triangular lattice one can find a simple analytical equation for the lowest 
instability frequency ui^kjj): 



(46) 



and for the current at the instability point 



Jz = 3 J 



V 



2wa ( v, 



+ ~ tI T- ) (^V 71 + ^i w A - VabUA) +ly 



The triangular lattice is stable at lo e < u>A(fcff)> which means that it always becomes unstable before reaching the 
resonance frequency lu p (it) = kn/2. The instability point approaches the resonance with decreasing dissipation. In the 
limit of high in-plane dissipation u a f, 3> 2/fcjj the instability point lu^ approaches the universal dissipation-independent 
value u>a — > kjiflypi = uj p (n)/\/2 

The stability boundary of the rectangular lattice uju{ku) is given by 



l 2 k 2 H = -WoV c (vabUn + + ^Ib^u 



In contrast to the triangular lattice, the rectangular lattice is stable at u E > u)n{kji) and remains stable at the 
resonance uj e = Iku- 
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C. Short-wave length instability 



To investigate stability of the lattice with respect to alternating deformations we substitute <p n — nn + u(— 1)" into 
Eq. (p5|). In the linear approximation with respect to the modulation u, it obeys the following simple equation 

du 

v 7r —+a 7r u = 0, (47) 

where 

a-* = a v (tt, k H , oje) = ~ Re [Q(k) - Q(tt - n)] (48) 
= v^(k, k H , uj e ) = v c + - Im [T (k - tt, k)] (49) 
are the "7r-stiffness" and "7r-viscosity" constants. The stability condition is given by an inequality 

o-kIv-k > (50) 

The "7r-stiffness" constant has a symmetry property ^(k, kH, we) = —a n (TT — k, for, u;.e). In particular, 
a^(Tt/2, ku, Wb) = 0, which means that the line k — tt/2 represents the stability boundary for alternating defor- 
mations at all lattice velocities. This means that no stability region can cross the line k = tt/2 and it is impossible to 
evolve continuously from the static triangular lattice to the rapidly moving rectangular lattice without intersecting the 
instability boundary. At small velocities the lattices with k < tt/2 are stable. Above a certain velocity situation is re- 
versed, the lattices with k > tt/2 become stable. This reversal point is determined by the condition dRe[G(K)]/dn = 
at k — tt/2. Any stability boundary in the region < k < tt/2 due to the sign change of o^(k) has the symmetric 
counterpart in the region tt/2 < k < tt. 



D. Stability phase diagrams 

In general, the lattice can become unstable at an arbitrary wave vector < q < tt. To find the stability regions 
of the moving lattice, we numerically scanned the real part of the decay rate a(q) (^) throughout the u>e — K phase 
diagram and find boundaries at which Re[a(<7)] changes sign at least for one q. The obtained stability phase diagram 
for the representative parameters v c = 0.002, = 0.1, and kn — 8 is shown in Fig. ||. For these parameters we found 
three stability regions at not very high lattice velocities (Josephson frequencies): the first low- velocity region is located 
below the resonance line and at tt/2 < k < tt, the second one is located above the resonance line and at tt/2 < k < tt, 
and the third high-velocity region is located along the resonance line at high velocities with n approaching with 
increase of cue- At the boundary of the first region the lattice experiences the long-wave instability for k > 2.04. At 
smaller k instability occurs at finite wave vector q = qi and the instability wave vector continuously grows with 
decrease of k. We find that this behavior occurs due to the sign change of the quartic term in the q-expansion of 
a{q) and consider in detail this transition in Appendix A. The instability at a finite q always means that the unstable 
mode has finite frequency tOi. Close to the instability point the system is expected to generate oscillations with this 
frequency. 

In the velocity range kn «C loe ;C knl there is only one stability region bounded by two lines above and below the 
resonance line ui p (k) fc# / \/2(l — cosk), ki(uje) < k < k u (^e), see Fig. |[ Both lines correspond to the long-wave 
instability. In a wide range of frequencies these lines are described by simple analytical formulae, which we derive in 
Appendix B: 

Ki (cje) = at \/v ab <uj e < (kHl) V3 /(Qvab) 1/3 (51) 

VVLJE 



Ku{ue) = k H J \J^~> B,tl/u ab < UJ E < v c 1/3 



(52) 



As one can see from Fig. ^, these asymptotics agree very well with the numerically calculated stability boundaries. 

Fig. ^ shows the evolution of the stability diagram with increase of magnetic field (parameter k}j, see Table 1). 
Salient features of this evolution are (i) the stability region above the resonance line shrinks with increase of field and 
vanishes at ku > 10, (ii) the high velocity region expands to higher k with increase of field, at ku ~ 14 additional a 
small stability island appears above the line k = tt/2, and at kn ~ 16 this island merges with the the low velocity 
stability region. This roughly corresponds to the field yj a a b / CcT^o / 2tt^s 2 . 
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VI. BOUNDARY CONDITIONS AND DYNAMIC STRUCTURE SELECTION 



A. General case of sharp boundary 



Stability analysis of the previous Section does not show what structure is actually realized at given velocity. For 
the static case the structure is selected by the minimum energy condition. Such condition is absent in the dynamic 
case. In this case a particular structure can be selected by the boundary conditions. Influence of the boundary on 
the moving structure is in turn determined by the interactions of electromagnetic waves with it. 

Consider a semiinfmite stack of junctions with n — 1,2... separated by a sharp boundary from a medium with 
arbitrary electromagnetic properties. We will obtain boundary conditions for such system in the case of a steady-state. 
To derive equations for the phase shifts <p n for such system we have to find solution of the linear equations without 
Josephson coupling taking into account boundary conditions. For the plasma wave with given frequency u> and wave 
vector along the layers k the oscillating phases (6 n ) and magnetic fields (h n ) in the junctions can be represented as 



'n ■ 



h n oc exp(— iq+n) + Bexp(iq + n), at n > 1 (53) 



where q + = q + (k,uj) is the complex wave vector given by Eq. (JL9J). Properties of the boundary are completely 
characterized by the complex amplitude of reflected wave B = B(k,uj), which has to be found by matching the 
solution ( |53] ) with electromagnetic oscillations in the medium at z < 0. In the case of weak dissipation at uje > fcff/2, 
the wave number q + has only small imaginary part and Eq. (^) describes the usual case of reflection of a propaga ting 
wave. On the other hand, at we < the wave number q+ is an almost pure imaginary number and Eq. (|53| ) 

describes reflection of a decaying wave, which is rarely considered in standard electrodynamics. Note that, in general, 
B(k,uj) can be a complex number with an arbitrary absolute value. Only in the simplest case of vanishing dissipation 
and propagating wave (Im(q+) = 0) the amplitude B(k,u) determines a conventional reflection coefficient, R{k,uo) — 
\B{k,uj)\ , and has property |B(fc,o;)| < 1. In the continuum limit, q+ <C 1, B(k,u>) is given by Fresnel formula (see, 
e.g., Ref. H) 

B(k,co)= £q+/ ;- £a f\ qt (54) 
where £ a fc(w) = e a m + 47 ^ ab — A /^ 2 is the in-plane dielectric constant of the superconductor, e is the dielectric 



constant of the medium at z < 0, and q t = y ^ k 2 is the wave vector of the transmitted wave. 

To derive equati ons f or the phase shifts valid close to the boundary one has to find the phase response function 



G(n,m) from Eqs. (16a, b) with appropriate boundary conditions. Due to relation ( |53| ) the functions G{n,m) should 
behave as 

G(n, to) oc exp (— iq+(n — ni)) + B exp(jg + (n + to)), at 1 < n < m. 



Solution of Eqs. ( [L6a| ,b) satisfying this condition can be constructed as 

G(n,m) = G(n- to) + BG(n + m), (55) 

where the second term describes the surface contribution and vanishes at tl , to — ► 00 . Equations for the steady-state 
phase shifts <j> n , valid near the boundary, have the following form 

_^ 00 

- 22 I m [G{n, to) exp (-i (4> n - 4> m ))] = ij (56) 

m— 1 

In general, a simple anzats <f> n — nn does not satisfy the steady-state equations (|56| ) near the surface. A general 
solution has the form (f> n = an + u n , where u n is the surface deformation, u n — > at n — > 00. Equations for u n can 
be represented in the form 

1 v 

- ^ Im [G(n, to) exp (— i«(n - to)) (exp (— i (u n - u m )) - 1)] = i s (n, n) (57) 

m— 1 

where 
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i s (n, n) — ij — — Im [G(n, to) exp (— «/t(n — m))] (58) 

m— 1 

is the excess Josephson current near the surface, i s (K, n) — > at n — > oo. The system ( |57j ) is degenerate because 
it contains only the differences u n — u m . As a consequence, the solution for u n exists only for certain values of k, 
i.e., bulk structure is selected by the boundary conditions. This mechanism of structure selection is realized in many 
dynamical systems (for a general review see Ref. |36j ). 

The condition for structure selection can be formulated in the explicit form in the case of small surface deformations 
u n <C 1, which is always realized at small velocities. In this case Eqs. ([37]) are reduced to a linear system 

^ oo 

— Re [G(n, to) exp (— iK,(n — m))] (u n — u m ) = i s (n, n). (59) 

TCI — 1 

The system (^9|) is degenerate and its solution exists only for special values of k.0 To formulate condition for the 
existence of the solution one has to find a solution ip n of the adjoint homogeneous equation 

oo 

(Re [G(n, m) exp {—in{n — to))] ip n — Re [G(n, to) exp in(n — to)] if> m ) = 0. (60) 

m=l 

Eq. ( |59| ) has a nontrivial solution only if its right-hand side is orthogonal to 'ipn, he., 

oo 

^^„i s (k, n) = 0. (61) 



Eqs. (|60| ) and (|61j) determine the lattice wave vector « at small velocities. 

For finite system, n — 1, 2, . . . , N, with identical boundaries the configuration is typically symmetric with respect 
to the midpoint n = N/2, i.e., 

ku, at 1< n < N/2 
k(N - n), at 1 < N - n < N/2 

In general, these solutions do not match at n = N/2, which means that they should be separated by a strongly 
perturbed region (shock). Further numerical simulations confirm such structure of steady-states in finite systems. 

Large class of boundaries is well described by the ideal reflection B = — 1. This includes a boundary with an 
insulator or free space (see Appendix C). In this case the response function G(n,m) in Eq. (56) acquires the form 

G(n,m) = G(n-m)-G(n + m). (62) 

We also calculate in Appendix D the reflection amplitude B{kf{,^E) for the more complicated but practically inter- 
esting case of the boundary between the static and moving Josephson lattices. Below we investigate in detail the 
evolution of the lattice structure for the ideally reflecting boundary. 



B. Lattice structure at small velocities 



At small velocities the lattice structure is close to the static triangular configuration and the phase shifts can be 
represented as 



0n( r ) = + u(t, n) 



(63) 



with \u(t, n + 1) — u(t, n)\ <C 1. Firstly, we find the lattice wave vector at small velocities from Eqs. (58), (|60|), and 
(|6l|). At lue — ► the function G(n) is real and it is only nonzero at n = — 1, 0, 1 with 



/•2 



G(0) «-_ 2 + - ; G(±l)«- r 



1 



/•2 



In this case the solution if) n of Eq. (|6(]) reduces to a constant and the condition ( |6l| ) yields 
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} J n) = 0. 



(64) 



For the ideally reflecting boundary, in a linear approximation with respect to u>e and x = 7r — k, the only nonzero 
components of i s {n, n) at n = 1 and 2 are: 



i a («, 1) 



2fc ff 



2^(3 + 1/Z 2 ) ^ w £ 

fab H 79 I 7772": 2 ) 



k 2 



2k 



H 



2k H 



Using these expansions we obtain from (p3) the following relation 



X 



7 + 2/1 2 

Vab H 72 | W E . 



(65) 



This relation determines the deformation of the lattice at small velocities for the ideally reflecting boundary. 

To obtain a general dynamic equation for weak lattice deformations, we substitute expansion ( |65| ) into Eq.(^5|), 
replace the discrete variable n by a continuous variable z, and perform a gradient expansion with respect to du/dz. 
This yields the celebrated Burgers equation for u(t, z): 



du 
'dr 



a 2 - 



d 2 u 



du 



Si, 



with 



«2 



6 2 = 



--Re 

4 

— Im 

4 



d 2 G(n;0) 
dn 2 
d 2 G(n; lu e ) 



2k 2 H 



v = v c H — Im 
2 



ck 2 
dG(ir, we) 



8u c \ uj e 



2k H 



dujE 



1 rr ' "■ 



lj e =0 V i? / n H 

Au c \ 2uj e 



k 2 ' 



Consider a steady-state configuration, du/ dr = 0. For a finite system with N junctions and identical boundaries the 
deformation is symmetric with respect to the center z = N/2. The regions of regular lattices, located at z < N/2 and 
z > N/2, are separated by the transitional region near the center where the deformation obeys the equation 



fu 



with 



a = 0,2/ 



lb, » \(v ah + 



UJE 



X = tt — k, and condition du/dz — > ±x at z — N/2 — > ±00. This equation has 



the analytical solution (see, e.g., Ref. pq) 



u = a In I cosh — { z 



a \ 2 

Using the result (|65|), we obtain that the width of the transition region L t — a/\ is given by 

1 



U = 



Vab+^) Ub+^W 



This length diverges at lue — > as uj^ 2 . This means that for a finite system there is a typical crossover frequency lun, 
which is determined by the condition L t (tON) — N. At lue < ujn the deformation field has a parabolic shape, typical 
for strained static systems. At uje > the system is split into the two homogeneously tilted lattices separated by 
the shock. 
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C. Numerical exploration of steady-states and structure evolution 



To find steady-state configurations at all velocities we solved Eq. (56) numerically for the case of ideally reflecting 
boundary using the same representative parameters (y c = 0.002, v a \, = 0.1). The evolution of dependencies k(oje) 
obtained from these solutions for kjj = 8 is shown in Fig. || together with the stability regions. We found that at small 
velocities (Josephson frequency u>e) the lattice structure smoothly evolves with increase of velocity and the lattice 
wave vector k decreases from 7r at zero velocity towards tt/2. The lattice structures in this region are shown in the 
left columns of Fig. 0. At certain velocity the lattice crosses the instability boundary. This velocity corresponds to 
the endpoint of the first flux-flow branch. It is close but not identical to the minimum velocity, ujE/kn — c m i n , of the 
plasma wave (in reduced units c m i n = 1/2 and in real units c m i n = cs / (2X a i,y/e^) w 3 • 10 7 cm/sec). At higher velocities 
the periodic lattice with the phase shift smaller than n/2 is restored. The structure continues to evolve smoothly 
towards a rectangular configuration with increasing velocity (see right column at Fig. ^) . 

Fig. |?j shows the evolution of the current-voltage dependence with increase of the magnetic field. At not very high 
magnetic field (kn < 16) the current-voltage dependencies have two stable branches, corresponding to the moving 
regular lattices. The low-velocity branch corresponds to the structure close to a triangular lattice. It terminates at the 
critical velocity. At small magnetic fields one can see a pronounced current enhancement prior the instability point. 
The high-velocity branch corresponds to the structure close to a rectangular lattice. The resistivity at this branch is 
close to the c-axis quasiparticle resistivity. In a wide range of magnetic fields, the stable branches are separated by 
the broad unstable region where a homogeneous regular lattice state can not exist. In this regime there is a broad 
range of currents, within which two lattice solutions exist. 

To characterize intensity of the electromagnetic wave generated by the moving lattice we plot in Fig. || the electric 
field dependencies of the Poynting vector along the layers (^) at different magnetic fields. The z component of 
the Poynting vector is always smaller than the ab component, and its sign corresponds to the energy flux from the 
surface to the bulk. One can see that the intensity of the electromagnetic wave rapidly increases with the electric 
field reaching a maximum at the instability point. This enhancement of the intensity of the traveling wave indicates 
that the increase of the current near the instability point is caused by the pumping of energy from a dc source into 
this wave. This maximum becomes smaller at higher magnetic fields. At higher electric fields the wave intensity 
decreases with field. Note however that one can expect another resonance at very high lattice velocities, where the 
velocity reaches the maximum velocity of electromagnetic wave, c max = c/e c (see Section [V). We do not consider 
this resonance here. 

To clarify the role of specific boundary on the forming the lattice structure and the current-voltage dependence, 
we also made several calculations for the more realistic case of the boundary with the static lattice. The amplitude 
of the reflected electromagnetic wave B(k,u>) for this case is calculated in Appendix D. Fig. || shows comparison 
of the structure evolution and the I-V dependence for the representative parameters v c = 0.002, v a \, = 0.1, and 
ktf = 6 with the case of the ideally reflecting boundary. Both cases show an overall similar behavior. The only 
qualitative difference exists at very small loe, where the the frequency falls within the range of the acoustic branch in 
the oscillation spectrum of the static Josephson lattice, < uj < \Z2/kn (see Appendix D). In this range the lattice 
wave vector exceeds it and the energy flux is directed from the bulk to the surface. However, this anomaly is almost 
invisible in the I-V dependence. It only becomes noticeable at smaller fields or weaker dissipation. In the case of 
the boundary with the static lattice the structure is closer to a triangular lattice at the first instability point. As a 
consequence, it has higher current at this point. This difference decreases at higher fields. 



VII. MULTIPLE BRANCHES DUE TO DYNAMIC PHASE SEPARATION 



We found that for not very high magnetic fields two stable states moving with different velocities may coexist 
within a certain range of applied current. The dominating in-plane dissipation strongly facilitates such coexistence, 
because in this case the lattice velocity is very sensitive to the lattice structure. Within the coexistence region one 
can expect a family of intermediate states, in which the system is split into the two regions moving with different 
velocities separated by a phase defect (dynamic phase separation, see Fig. JlC|)), In a continuous system coexistence 
of the states moving with different velocity at the fixed average velocity is possible only at one driving force, at which 
these states are in the dynamic equilibrium. Out of equilibrium the boundary will move and eliminate one of the 
states. However, for our discrete system the boundary separating different lattices is pinned by the discrete structure, 
and coexistence is possible in a wide range of driving forces (transport currents). 

If the fraction f s of the system is in the slowly moving state with velocity v s and the remaining fraction 1 — f s is 
in the rapidly moving state with velocity Vf than the average velocity v av , which determines the observable electric 
field, is simply given by 
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Vav = f S V s + (1 - fs)Vf 

To investigate the structure of the boundary region, consider an infinite stack, in which the rapidly moving state 
occupies the region n > and the slowly moving state occupies the region n < 0. The condition for coexistence of 
such states is given by 

v c ujf + ij(ujf,Kf) = v c uj s + ij(uj s ,K s ) 

where ujf = kjjVf (uj s ~ knVs) is the Josephson frequency for the fast (slow) state. In the zero order approximation 
with respect to the Josephson coupling the phases are given by 

9 n = LUfT + kiju + Kfn + u n , at n > 
= lu s t + knu + n s n + v n , at n < 



To obtain equations for the boundary deformations u n and v n we have to find oscillating phases induced by the 
Josephson currents. The phases and fields, oscillating with the frequency u>f, are determined by equations 
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ab dT \ 
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where Q(n) is the step function. Solution for the oscillating phase is given by 



= Im 



E 

m-0 



G f{n — m) exp (i (w/r + knu + KfTi + u n )) 



with Gf(n — to) = G(n — m; Kf, ojf). Equation for u n is obtained from condition (sin# 

Im[G/(n — to) exp (i (Kf(m — n) + u m — u n ))] 
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m— 
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A similar equation can be derived for the boundary deformations of the slow state v n . Comparing this equation with 
Eqs. ( p5| ) and (|56|), one can see that the deformations at the boundary between different moving states coincide with 
the deformations near the nonreflecting boundary (B = 0). 

The phase-separated states give the most natural interpretation of the multiple I-V branches observed by JJecht- 
fischer et al. who studied transport properties of Josephson lattice in BSCCO mesas at high magnetic fields.0 Note 
that these branched should not be mixed with the multiple branches at zero magnetic field, which appear due to 
the switching of the separate junctions into the resistive state. This interpretation can be verified by measuring 
the spectrum of the microwave irradiation emitting by the stack. Instead of a single peak located at the Josephson 
frequency corresponding to the average voltage, the spectrum of irradiation should contain two peaks: at the "slow" 
frequency uj s with the weight f s and at the "fast" frequency to / with the weight 1 — f s . 



VIII. CONCLUSIONS 



In summary, we performed a detailed investigation of the Josephson lattice dynamics in layered superconductors 
at high fields. Our main results can be summarized as follows: 

• Interaction between Josephson planar arrays in layers is mediated by exited electromagnetic oscillations. The 
dynamics of these arrays is described by nonlocal and nonlinear equations. 

• For the coherent steady-states, the excess Josephson current and the Poynting vector of the generated electro- 
magnetic wave have resonance dependence on the Josephson frequency (electric field) . The resonance frequency 
is the plasma frequency at the wave vector selected by lattice structure. In the case of the dominating in-plane 
dissipation, typical for high-temperature superconductors, the damping parameter also strongly depends on the 
lattice structure. 
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• We investigated stability of the coherent steady-states and found the two major stability regions in the plane 
lattice structure-lattice velocity: the low-velocity region and the high-velocity region. Exact topology of the 
stability phase diagrams depends on dissipation and magnetic field. 

• Lattice structure at given velocity is selected by boundary. Boundary conditions are determined by the reflection 
properties of electromagnetic waves, generated by the moving lattice. 

• We investigated structure evolution with increase of the lattice velocity. In a wide range of fields there are two 
stable branches, the low-velocity branch and the high-velocity branch. At low velocities the lattice structure 
is close to a triangular lattice. This low-velocity branch terminates due to instability at the critical velocity 
near the minimum velocity of electromagnetic wave. At high velocities structure of the lattice approaches a 
rectangular configuration. 

• Experimentally observed multiple branches in I-V dependencies arc interpreted as the phase separated states, 
in which system is split into the slowly and rapidly moving regions. 
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Appendix A: Continuous transition from the long-wave instability to the finite-q instability 

Numerically investigating stability of steady-state we found that the nature of instability of the first flux flow branch 
depends on the wave vector k, at which the lattice becomes unstable. If k is larger than the critical value Kt than 
instability occurs at q = 0, i.e., d 2 ai/dq 2 changes sign. At k < Kt the lattice becomes unstable at the finite wave 
vector qi, i.e., ax(qi) = at the instability point. The transition is continuous, qi smoothly grows with decrease of k 
starting from qi = at k = Kt (see Fig. [Ly). This behavior can be explained as follows. At small q the decay rate 
a.\ (q) can be expanded with respect to q 

a(q) w -a 2 q 2 - - ^a 6 <? 6 . 

The coefficients a 2n in this expansion depend on k and We- For stable lattices a 2 > 0. The coefficient a 2 becomes 
negative when the frequency ui E exceeds the long- wave instability boundary lo 2 {k). Near w 2 (k) one can use linear 
expansion for a 2 

a 2 = ci 20 (lj 2 (k) - we)- 

Lets also define the frequency w±(k), at which the coefficient changes sign. Near w^{k) we have 

a 4 = aio^fft) - We)- 

If wi(k) > lu 2 (k) then the instability occurs at q = 0. For the first flux flow branch this is the case for k > K t . 
Otherwise the lattice become unstable at a finite wave vector qi. The transition between these regimes takes place at 
the lattice wave vector K t defined by 

W 2 (K t ) = 0J4,(K t ), 

i.e., both a 2 and 04 vanish at k = Kt- Lets consider in detail the behavior in the region where k only slightly smaller 
than Kt and the difference w 2 (Kt) — w^Kt) is small. In this region we can use linear expansions for both a 2 and 04. 
We define 

5 = lo e - W4,(k), 

5 2 = w 2 (k) - cj 4 (k) = (ui' 2 - w' 4 ) (k - K t ) > 0, 

and represent a(q) in the form 
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a(q) = -a 2a (S 2 - 8)q 2 + -a^Sq 4 



r *q 



With increase of 5 the dependence a(q) becomes nonmonotonic at (a^S) > 4a2o&6(<52 — 5). Extremum values of q 
are determined by solutions of quadratic with respect to q 2 equation 



a 6 q 



a i0 5q 2 + a 20 {6 2 - 5) = 



(66) 



and are given by 



4 




where q_ gives the minimum of a(q) and qi_ gives its maximum. The lattice becomes unstable when cx(q+) = 0. This 
equation together with the extremum condition ( |66| ) determines the instability point 8 — Si and the wave vector of 
the unstable mode g. 



8 06020 




Close to the transition 8 2 <C 06^20/040 we have 



16 a 6 a 20 
3a4o8 2 
4a 6 



40 S 2 
°2 i 



This means that the wave vector of the unstable mode grows as g, oc n t — k near the transition. The calculated 
dependence qi (k) agrees very well with this predictions (see the lower-left inset in Fig. . 



Appendix B: Stability boundaries in the intermediate frequency range 

In this appendix we consider the stability region in the intermediate velocity range, <C lje <C k^l, where it is 
possible to derive simple analytical expressions for the stability boundaries. There is only one stability region in this 
frequency range located at k <C 1 and limited by two stability boundaries, ki(u>e) and k u {u>e), located below and 
above the resonance line. Both bo unda ries ki(u>e) and k, u (uje) correspond to the long-wave instabilities and we will 



use general relations from Section V B . Assuming inequalities l/l 2 , h> a b^E/l 2 < k 2 « 1, v c k 2 -C v a b and introducing 
scaling variables 



z = k UJ E /k H , y = v ah io E , 



we present the coefficients from Eqs. (42a-d) in the form 



a 2 



k 2 



2k H 



Im 



Re 



v = v c k H H — 

2(j}Ekn 



C = —ujEkntC Re 



(z-1- iy)' 

(1+jy) (3z + l + iy) 
(z-1 - iyf 

—2z + iy 



Im 



(z- 1 -iyY 
2 + iy 



(z-1 - iyf 



(67a) 
(67b) 
(67c) 
(67d) 
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In general, the stability boundaries are determined by Eq. However, at high enough Josephson frequency term 

aiQ/u can be neglected. We focus here on the regime i/ > 1. Consider the lower boundary ki(we) first. As can be 
checked from the result for this boundary the term aiQ/v is smaller than a?, at least by the factor 1/y. Assuming 
inequality i/>zwe can approximate 0,2 in Eq. ( |67b| ) by the first expansion term with respect to 1/y, di ~ — ■ 
Therefore the lower stability boundary is simply given by z = 1/6, consistent with the assumed condition i/>z. This 
is equivalent to 



ki (w E ) 



k 



H 



V6t 



we 



(68) 



This expression for the lower boundary is valid up to we < (knl) 2 ^ 3 /(6t / a fc) 1 ^ 3 ■ It is interesting to note that this 
boundary does not depend on the dissipation parameters. 

Let us consider now the upper boundary k u (we)- The term a\(,/v in Eq. (Q) can again be neglected if k 2 <C 
UJ E^H v cV a b- In this case the upper boundary is again determined by the equation ai = 0. At y 3> 1 this equation has 
asymptotic solution in the form z = ay 3> 1 and in the main order with respect to 1 jy the constant a is determined 
by equation 



Re 



i (3a + i) 



= 0, 



which gives a 



(a - if 

'3/5. This corresponds to the stability boundary 



k u (we) = k H \ \ ~ — 
\\5w E 



(69) 



Substituting this expression into the inequality n 2 -C w%kjji/ c v a b (condition to neglect the term aiQ/v in Eq. (p4[)), 
we see that it is equivalent to the condition we S> v c which gives the lower limit of applicability of Eq. (|69|) . 



Appendix C: Boundary with an insulator 

Consider a semiinfinite stack of junctions with n = 1,2,... bounded by an insulator at semispace z < (e.g., 
by free space). We assume that the transport current is fed to the stack through the first layer. In this case the 
oscillating electric and magnetic fields induced by the moving Josephson lattice should match with the electromagnetic 
oscillations in the insulator. 

The oscillating electric and magnetic fields decay into the insulator as 

E(r, t) — Eo exp (iknx + qz + iwst) , 
H(r, t) = Hq cxp (iknx + qz + iwst) , 



where E = (^0,0,^0), H = (0, H Q ,0), q — \/k 2 H ~ ew E /c 2 , and e is the dielectric constant of the insulator. 
Maxwell's equation connects E x q and Hq as 



-qH Q 



iwe£ 



E x o- 



Continuity of the parallel component of the electric field gives the relation 



$ iw 

2-KC 



and Eq. (g) gives 



-IWE 



$0 

2ttX 2 . 

ab 



pi 



From the last three equations we obtain the boundary condition 

e(X ab w E /c) 2 



k 2 H ew 2 E /c 2 H Q = - 



Hi - H 



Air\ 2 <j ab iw E /c< 



{Hi - H a ) 



(70) 
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In general case one has to use the complex dielectric constant e of the medium at the frequency u>e and wave vector 
k[{ . In reduced variables this equation can be rewritten as 



k 2 H -—^cj%ho = — — ^ (hi -h ) (71) 

ej 2 a e c 7 1 + iv ah w E 



From this relation we obtain the amplitude of reflected wave 



„ J k H ~ + ^T+fel (1 - exp(-. g+ )) 

In the case e ~ 1, fc# ~ 1/7*1 an d ~ cj p we obtain the estimate /iq ~ h\j^e c -C /ii. Therefore the condition ( |7l| ) 
with high accuracy may be replaced by the much simpler condition h — 0, which corresponds to the ideal reflection 
case, B = — 1, 



Appendix D: Reflection amplitude for boundary between moving and static lattice 

Small mesas are typically used in experiment to enhance transport current density. In this case the transport current 
flowing in the bulk part of the sample is not sufficient to drive the Josephson lattice there and the moving Josephson 
lattice in mesa neighbors with the static lattice in the bulk (see inset in Fig. |l2]). In this appendix we consider 
electromagnetic properties of such boundary. We assume that the Josephson lattic e is static at n < and moves with 



constant velocity at n > 0. According to the general recipe outlined in Sec. VI A we have to match electromagnetic 
oscillations in the neighboring medium neglecting Josephson currents in the region of the moving lattice and to find 
the amplitude of the reflected wave. We consider reflection of the wave 9 n oc exp (i (knu + lot — iq + n)) coming 
from above with q + given by Eq. (|l9|). This incident wave excites phase and field oscillations of the static lattice. 
Such oscillations have been theoretically studied in Ref. |[ At high in-plane field these oscillations are described by 
approximate equations 

dh 

(uj 2 — w c lo) 9 n — (2C sin 2 (fcyu + mi) + cos (fc#u + wn)j 9 n H — = 0. (73a) 

ou 

7 2 1 \ , de n t . f d6 n h n 



Vi- T2 )h n + ^ + iu ab u^-f) =0, (73b) 

where C = ^cos 9n\u)\ w 4 ^/' and 9n \u) w k^u + nn — 4H fc V sin(fc//u + 7rn) is the static phase difference. The 
cosine term in the first equation couples oscillations with the in-plane wav e vect or kn to the homogeneous mode 



(at kn ^> 1 coupling to the higher harmonics can be neglected). Eqs. (73a) and ( |73q ) have two types of solutions 
for the oscillations with the wave vector k^, with 9 n oc cos(k^u) and 9 n oc sin (kuu), which we refer to as the 
even and odd modes. The incident wave excites both modes but only the even mode is coupled to the homogeneous 

oscillations of the lattice. For this mode we look for solution in the form 9 n (u) = (^9 cos (knu) + (—l) n (fj exp(— zx e n), 
h n (u) = h sin (fcjjti) exp(— ixe n ) and derive equations 

(oj 2 - iv c uj) 9 + k H h = 6, (74a) 
(u 2 -C -w c lo)9 = ®-, (74b) 
2(1 - cos Xe ) + l + lVabUJ \ h + k H (l + iv ah uj) 9 = 0, (74c) 



I 2 

which determine the eigenvalue Xe for the given frequency lo and in-plane wave vector kn as 

(u 2 - iv c uj -C) (1 + iujv ab ) k 2 H l + iujv ab 
2(l-cosx e ) = — ; — : — r . (75) 

(UJ Z — IV C UJ) (LU Z — IV C UJ — L) — 2 ' 

Being inverted, this equation gives the frequencies of plasma modes at fixed c-axis wave vector x and their decay. 
Note that at u> = we obtain the solution corresponding to a homogeneous shift of the lattice, \ — 7r - For weak 
dissipation this equation has two other important solutions: (i) w ~ VC ~ \/2/kH and Xe near zero (the endpoint 
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of the acoustical branch of the oscillation spectrum) (ii) u> « kjj/2 and \e ~ tt (the frequency of the homogeneous 
plasma mode). 

For the odd mode we look for solution in the form 9 n = 9 sin (fcjjtt) exp(— ixc.n); h n — h cos (kjju) exp(— zx Q n) and 
derive 

(w^g — i^cCij — C) 9 n — ktf h = 



( 2(1 - cosxo) + 1 + l !f Vab )h+(l + iLuis ab ) k H 9 = 



which gives 



2(1 - cos Xo ) = (1 + iuv ab ) ( 2 _ k " _ r -^) (76) 



Because of cos(fc#u) modulation in the static lattice region the reflected wave contains a backscattered component, 
with the wave vector — kn, so that the full solution for the oscillating magnetic field has the form 

h = exp(i/cifit — iq+n) + B exp(iknu + iq+n) + Bb exp(— ikuu + iq+n), at n > 
= A e cos kuu exp(— ixen) + i*4 sin/c#uexp(— iXori), at n < 

Introducing the reflection amplitudes in the even and odd channels, B e — B + Bb and B Q = B — Bb, we obtain by 
matching the solutions at n — 0, 1 the following equations 

1 + B e = A e 
exp(-iq + ) + B e exp(ig + ) = A e exp(-zx e ) 

and similar equations for B and A . These equations can be easily solved giving 

exp(-ixe) - exp(-i<7+) 



B e = - 



exp(-zxe) - exp(ig + ) 

2i sin(g_|_) 
exp(-ix e ) - exp(iq + ) 



and similar formulas for B and A Q . These formulas resemble the well-known Fresnel formulas of classical electrody- 
namics. The am plit ude of the reflected wave B = B{kn, oj), which determines the surface contribution to the response 
function in Eq. (p5|), is given by 

B= \( exp(-%) ~ exp{-iq + ) | exp(-zx Q ) - exp(-iq + ) \ 
2 V exp(-zx e ) - exp(zg+) exp(-?Xo) - exp(i<j + ) / 

Where, again, q + , Xe, and Xo are given by Eqs. jl9|), ( ^), and (76). Plot of the frequency dependence of the 
amplitude and phase of B(kn,uj), B(kn,uj) = exp(i0g), is shown in Fig. [l^. Typically \B\ is large in the small 
frequency range < uj 2 < C, where the incident decaying wave excites the propagating acoustic waves in the static 
lattice. At uj 2 > C \B\ rapidly goes to zero, because both media have identical spectrums of electromagnetic waves at 
high frequencies. Both the amplitude and phase of B(kn,oj) have anomalies at the typical frequencies u> — \f2jkn 
and u> = fcjj/2 of the oscillation spectrum of the static lattice. These anomalies become more pronounced at lower 
dissipation and smaller field. 
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Table 1. Meanings, definitions and practical formulas for the reduced parameters used throughout the paper. In 
practical formulas f p = uo p /2-k means plasma frequency, p c and p a b are the components of the quasiparticle resistivity 
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Notation 


Meaning 


jjcnnition^Ovjo ) 


rra-ctiCcu iorniLiia ^r}ov_yOvjj 


UJ E 


reduced Joscphson frequency 


L'KCShjz 


l c fP \ ItyVVI 

[s±!j z )[m\ j 




9 1 0— 3 f PTTv 


k H 


wave vector of Josephson lattice 


Z7TI1 
<t>„ 






c-axis dissipation parameter 


4ttct c 


1.8-10 a 




£r/9r- [Sl-Cml fnlGHzl 


Vab 


in-plane dissipation parameter 


47TCT ab A^ b W p 


0.79(A afc [ M m])^/ p [GHz] 


c 2 




I 


reduced London penetration depth 


Kb/ s 




h 


reduced local magnetic field 


2njXi b H 





(a) 
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♦ ♦ ♦ ♦ ♦ 
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(b) 



FIG. 1. Steady-states of moving Josephson lattice (squares mark positions where the interlayer phase difference is equal to 
7r + 2irm): (a) a periodic lattice, the structure is determined by the lattice wave vector k, which is connected with the shift 
between the lattices in the neighboring layers 8a by relation k = 2n8a/a (b)a double-periodic lattice 
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FIG. 2. The frequency dependence of the decay length L s for different values of the magnetic wave vector kH (see Table 
1) for the typical parameters v a b = 0.1 and v c = 0.002. This length determines the interaction range between the Josephson 
planar arrays in different layers. A finite stack, containing N junctions, shows bulk behavior if N 3> L 3 . 
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FIG. 3. Stability regions of the moving Josephson lattice in the plane ue-k calculated for the representative parameters 
u c = 0.002, j/ a {, = 0.1, and kn = 8. Grey regions correspond to unstable lattices. Sections of boundaries marked by black 
correspond to the long-wave instability. Sections marked by light grey correspond to the short-wave instability. Line starting 
at (0, 71") shows dependence of the lattice wave vector on the frequency uje (lattice velocity) selected by the ideally reflecting 
boundary. We also show the resonance line, corresponding to matching between the Josephson frequency and the frequency of 
the plasma wave at the wave vector k. 
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FIG. 4. The boundaries of stable region in the frequency range ku <SC we -C kul for the same parameter set as in Fig. ^. 
The line in the middle is the resonance line. Dotted and dashed lines show the asymptotics of the lower and upper boundaries 
given by Eqs. ( [il|) and . 
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FIG. 5. Evolution of the stability diagram with increase of the magnetic field (for 7 = 500 H ~ ku ■ 0.3T). 




FIG. 6. Lattice structures in the real space at different values of the Josephson frequency loe (lattice velocity). The structures 
are generated using numerically calculated phase shifts for the system with 51 layers and the same parameters as in Fig. [| 
(y c — 0.002, v a b — 0.1, and ku = 8). The left column shows the lattice structures in the low-velocity stability region, below 
the first instability point. The right column shows the lattice structures in the high-velocity stability region. 
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FIG. 7. The current-voltage characteristics at different magnetic fields, for ku = 6, 8, 12, and 16. The dependencies are 
obtained using numerically computed steady-states with parameters v c — 0.002, i/ a t = 0.1. Thick lines show the stable branches 
and thin lines show the unstable branches. The branches, corresponding to the double-periodic states, are marked by dashes. 
The plots are displaced vertically for clarity. 



0.2 




2 4 6 8 10 12 



FIG. 8. Electric field dependencies of the Poynting vector along the layers for electromagnetic wave generated by the moving 
lattice (see Eq. (|3^)). The scale P a b is defined by Eq. (^). We use the same notations for different branches as in the previous 
figure. 



2G 



0.08 
>0.06 
0.04 
0.02 



(b) 




- ideal reflection (stable) 

- ideal reflection (unstable) 

- bound, with static lattice 



10 



12 



CO 



FIG. 9. This figure illustrates influence of the boundary conditions on the structure evolution (a) and the current- volt age 
dependence (b). We compare the ideally reflecting boundary (B = —1) and the boundary with the static lattice (B(k,u) is 
calculated in Appendix D). The comparison is made for v c = 0.002, v a b = 0.1, and ku = 6. Both cases show an overall 
similar behavior. In the second case there is a region at small Josephson frequencies, < uj < \/2/kH, of the anomalous 
structure evolution, where the lattice wave vector exceeds n (this region is blown up in the inset). This region corresponds to 
the frequency range of the acoustic branch in the oscillation spectrum of the static Josephson lattice. The dependence k(uie) 
has a kink at the endpoint of the spectrum. However this point is almost invisible in the I-V dependence. The second case is 
also characterized by the stronger current enhancement near the instability point (we show only stable I-V branches for this 
case) . 
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FIG. 10. Multibranch structure of the current-voltage characteristic due to the dynamic phase separation. Two states 
corresponding of the slow lattice motion (velocity v s ) and the fast lattice motion (velocity Vf) coexist within the current range 
marked at the vertical axis. In this region the intermediate phase-separated states exist, in which the system is split into the 
rapidly and slowly moving regions. The intermediate branch corresponding to one of such states is shown by dotted line. 
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FIG. 11. Behavior of the wave vector q; (left axis) and frequency u>i — Im[a(qi)] (right axis) of the unstable mode. The 
upper-right inset shows the transition point in the k-uje diagram of Fig. ^| Lower-left inset shows plot qf (k) with linear fit 
below the transition. 
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FIG. 12. The frequency dependence of the absolute value (\B\) and phase (<f>s) of the amplitude of reflected wave for the 
boundary between moving and static lattices (Eq. (|t7|)). The ploted curves are computed using parameters v c = 0.002, u a t = 0.1 
for two values of ku, 4 and 8. The typical frequencies of the oscillation spectrum of the static lattice ( to = y/2/kH and ku /2) 
are marked for kn = 4. Inset sketches a small mesa on the top of bulk crystal. Concentration of the c-axis transport current 
inside the mesa forces the lattice move only in this region. 



28 



